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In qubits made from a weakly anharmonic oscillator the leading source of error at short gate times 
is leakage of population out of the two dimensional Hilbert space that forms the qubit. In this paper 
we develop a general scheme based on an adiabatic expansion to find pulse shapes that correct this 
type of error. We find a family of solutions that allows tailoring to what is practical to implement for 
a specific application. Our result contains and improves the previously developed DRAG technique 
[F. Motzoi, et. al, Phys. Rev. Lett. 103, 110501 (2009)] and allows a generalization to other 
non-linear oscillators with more than one leakage transition. 



I. INTRODUCTION 

The physical realization of quantum information pro- 
cessing in superconducting circuits [1-4] has enjoyed re- 
markable progress over the last decade. While initially 
decoherence limited single qubits to only a few coherent 
oscillations [5] , high precision, general quantum conrol is 
now possible over single- and fcw-qubit systems. This is 
evident by the demonstration of high-fidelity nonclassical 
states of two-qubit [6-9] and three-qubit [10, 11] systems, 
harmonic oscillators [12], and the demonstration of small 
quantum algorithms [8]. 

This success is partially due to our current understand- 
ing of sources of noise and the development of techniques 
and systems that are resilient to these noise sources. Ex- 
amples include the optimum working point [13] and the 
introduction of low-dispersion qubits like the transmon 
[14, 15] and the capacitively shunted flux qubit [16]. On 
the other hand, a promising route to success are qubits 
that contains only a minimal number of elements, such 
as the phase qubit [17-19]. What these systems have in 
common is a weakly anharmonic energy level structure, 
i.e., the states that are outside of the qubit subspace 
spanned by |0) and |1) are only separated from each other 
and the qubit subspace by energies only slightly different 
than the qubit frequency. 

Having a weakly anharmonic qubit poses a challenge 
in the implementation of quantum gates. It is known 
[20, 21] that the time evolution operator of a linearly 
driven harmonic oscillator is a combination of a coherent 
displacement operator tracking the classical trajectory 
of the driven oscillator and a global phase factor. This 
evolution encompasses all energy levels and cannot be re- 
duced to a single-qubit rotation. Spectroscopically, this 
can be understood as follows: a single qubit rotation is 
typically implemented by a pulse of radiation resonant 
with the qubit energy splitting. In a harmonic oscillator, 
all energy splittings are the same, so driving one transi- 
tion drives all others at the same time. A system starting 
initially in an energy eigenstate will quickly be driven into 



a superposition over many energy eigenstates. By this to- 
ken, it is crucial that a qubit is nonlinear [1, 3, 4, 13, 22], 
that is that the transition frequency of the qubit lev- 
els is different by an amount A from the transition fre- 
quencies to the non-qubit levels. Spectroscopically, we 
would expect that whenever the bandwidth of the pulse 
comes close to A, i.e., when its duration becomes short on 
the scale of 1/A, we expect significant leakage to higher 
states. Thus, it is a challenge to implement fast single- 
qubit gates in weakly anharmonic systems. The imple- 
mentation of faster gates is important as it allows more 
gates to be executed in a given coherence time, an im- 
portant step toward high-fidelity quantum logic. 

While superconducting qubits are the most well known 
example of qubits made from weakly anharmonic oscilla- 
tors, there are many other examples. In fact, it has been 
shown that no physical particle can be a true qubit [23]. 
Examples of leakage states include higher vibrational 
states in optical lattices [24], polarized spin states in the 
singlet-triplet qubits [25], and auxiliary states in ion traps 
[26] and Rydberg atoms [27-29]. 

In this paper, we outline a suite of strategies to im- 
plement single-qubit quantum gates in qubits singled 
out from the spectrum of an anharmonic oscillator. We 
develop an adiabatic expansion technique that leads to 
order-by-order constraint equations on a toggling-frame 
transformation and the control fields. The space of solu- 
tions contains the DRAG strategy proposed in Ref. [30] 
and re-analyzed in the presence of an oscillator bath in 
Ref. [31] as a special case, as well as simpler versions of 
DRAG that do not require time-dependent energy bias 
or phase ramping, similar to those implemented in Rcfs. 
[32] and [33]. Additionally, we derive optimal solutions 
to a given order by minimizing the errors manifested by 
the fields in the next higher order. 

The outline of the paper is as follows: in Sec. II we 
review qubits made from anharmonic oscillator and intro- 
duce their rotating-wave description. Sec. Ill describes 
Gaussian pulses and shows how they lead to both phase 
and population errors. In Sec. IV we extend the DRAG 
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scheme of Ref. [30] to find a wealth of different pulse good approximation because parity forbids all other tran- 

shapes that give improved performance over Gaussian sitions. The control Hamiltonian is 
shaping. In Sec. V we apply the generalized scheme to 

cases where there are more then one leakage level. Finally 'tz^ 

we conclude in Sec. VI. H ct (t) = £(t) ^j-i a j-i,ji ( 2 - 2 ) 



II. THE SYSTEM 

A. Lab frame Hamiltonian 

We consider a qubit formed by the two lowest levels 
(which we generalize in sec. V) of an anharmonic oscil- 
lator. These levels are separated in energy by hui, where 
ui is the transition frequency. The j th higher levels are 
different to Hjuj by tiAj 1 where Aj is known as the an- 
harmonicity. That is, the Hamiltonian for the nonlinear 
oscillator of dimension d is (H = 1) 



d-l 



(2.1) 



where EL,- = \j)(j\ is the projection operator onto the j th 
energy level. Without loss of generality we set Ai = 0. 
This is illustrated in Fig. 1 with the qubit levels being |0) 
and |1) (green). For many nonlinear oscillators the an- 
harmonicity takes the form Aj — A2Q — l)j/2, which we 
will call the standard nonlinear oscillator (SNO), essen- 
tially a Duffing oscillator within the rotating wave ap- 
proximation [22]. For the lowest few levels, supercon- 
ducting qubits of the transmon [15], phase qubit [19], 
and capacitively shunted flux qubit [16] types are well 
approximated by a SNO. Furthermore, motional states 
in optical lattices [24], collective modes of ion traps [26] 
and nanomechanical oscillators [34] are also described as 
SNOs. 
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FIG. 1. (color online) Energy level diagram of the system 
we are considering. The qubit is formed by the |0) and |1) 
(green) levels and we aim to have complete control in this 
subspace when leakage to the |2) and then |3) etc. is possible 
(red arrows). The dotted black lines indicate the positions 
the energy levels would be at if the system was a harmonic 
oscillator of frequency u. 



where £ (t) is the drive amplitude, erj fc = |j)(fc| + 
is one of the effective Pauli spin operators for levels j 
and k, and Xj is a dimensionless parameter that weighs 
the relative strength of driving the \j) — > \j + 1) tran- 
sition versus the |0) — > |1) transition. In our model we 
take Ao = 1 and leave the Aj's as input parameters. For 
a harmonic oscillator controlled via a dipole interaction 
with an external field Xj = ^/j; however, in appendix A 
we show that in cavity or circuit QED architectures Xj 
can differ substantially from this value. 

For the functional form of the drive £ (t) we will assume 
that \£ (t)\ <C lj (weak driving regime) and introduce en- 
velope shaping of the driving field at carrier frequency 
LOd- This leads to 

£(t) = n x (t) cos(ujdt + 4>o) + Qy(t) sm(ojdt + <j)o)- (2-3) 

As per convention, the two quadratures amplitudes £l x (t) 
and fly (t) can be amplitude modulated using a waveform 
generator and then mixed back together with the carrier 
to give this form of control field. Here, 4>q is the relative 
phase between the envelope and the carrier at the start 
of the operation. This phase is irrelevant if the rotating 
wave approximation can be made, as will be shown in 
the next section. 



B. Rotating frame Hamiltonian 

For quantum information processing it is highly suit- 
able to define operations with respect to the frame ro- 
tating at the driving frequency Ud- In this frame we 
have three independent controls: 5(t) — uj(t) — u>d (the 
qubit detuning), Vt x (t) and S\(£), which, projected to 
the qubit subspace, control application of the three Pauli 
spin operators Uo,ii cr o.i' an d CT oi' respectively For ex- 
ample, the identity operation is achieved by setting all 
the controls to zero. Note that here we have assumed 
that 5 is controlled by shifting the qubit frequency. This 
is not necessary and, as shown in appendix B, this can 
be achieved by a time dependent phase in Eq. (2.3). To 
move to the rotating frame we define the unitary 



d— 1 



R(t) = ex P [-ijudt] n, 



(2.4) 



We will assume that control in this system is due to 
some dipole-like interaction that only allows single pho- 
ton transitions. As for harmonic oscillators, this is a 



which determines the transformed Hamiltonian 



H R (t) = R\t)H{t)R{t) + itf(t)R{t), (2.5) 
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where H (t) = + H ct (t). Explicitly we have 



d-l 



p -iLu d t-icf)a i ( 



xEVi[ij)(j-ir f +t.c] 



(2.6) 



where h.c. stands for hermitian conjugate and Q(t) = 

n x (t) + in y (t). 

Assuming that ujd is larger then any other rate or fre- 
quency in this frame we can perform the rotating wave 
approximation (i.e. time average the fast rotating terms 
to zero) . For the SNO case this amounts to restricting the 
dimension d to be less then •v/zcJ/Aj specifically d = 7. 
After this approximation we can write the Hamiltonian 
as 



ff fl (t) = ^(i«5 + A j )n 

3=1 
d-l 



(2.7) 



|i)(i-l|+h.c 



Here we have included the 4>q into the energy states 
~ * e 1 ^ \j)), and we see that within the rotating wave 
approximation the relative phase between the envelope 
and the carrier at the start of the operation is irrelevant. 
Finally, Eq. (2.7) can be rewritten as 



d-l 

3=1 
d-l 



(2.8) 



where a v - k — — i\j){k\ + i\k){j\ for k > j. We see that if 
we can restrict the system to the lowest two levels then 
all rotations in the single qubit space can be achieved 
by independent controls; however, in general, this is not 
true. In Sec. Ill we show that the higher level transitions 
lead to a combination of a phase and leakage error [30, 
35]. This has been experimentally measured in Refs. [36] 
and [37]. 



III. GAUSSIAN SHAPING AND ERRORS 

Our goal is to implement gates contained within the 
qubit subspace. That is, we want to shape Q x (t), & y (t) 
and S(t) in Eq. (2.8) so that 



Uideai = Texp 



-l I H H (t)dt 
o 



= e^£/ qb ©[/ rcst , (3.1) 



where t g is the gate time, T is the time ordering operator, 
f/qb is a unitary that acts only in the qubit subspace, [/ rcs t 
acts only outside of the qubit space, and 4> describes a 
relative phase. Therefore, J7 rcs t as well as the phase <f> are 
completely irrelevant for operations in the Hilbert space 
formed by the qubit. 

To demonstrate the typical set of errors we choose 
U q b — 0"oi, the NOT gate. For a leakage-free qubit 
this would be implemented by simply setting S(t) = 0, 
&x{t) = ^g(^)j an d Qy(t) — 0, with the only requirement 
that that J^ a Qcitydt = tt. To reduce the leakage to the 
third level, the standard result prior to Ref. [30] was to 
use Gaussian modulation of the envelope [35, 38]. In this 
case Og (t) takes the form 



Q G (t)=A 



exp 



(t-t g /2) 2 



exp 



~8<7 2 



V2na 2 erf[t g /VS<T} - t g exp[-t 2 J8a 2 } 



, (3.2) 



where t £ [0,t g ], a is the standard deviation and A is 
chosen such that the correct amount of rotation is imple- 
mented (e.g. A = 7r for a NOT). This functional form is 
chosen to enforce that the pulse start and end at zero. 
In the limit that t g — > oo we recover the standard Gaus- 
sian function. The motivation for Gaussian shaping is 
that the small and strictly limited frequency bandwidth 
(1/er) ensures little excitation at the leakage transition 
frequency. For short pulses however, there is still signif- 
icant spectral weight at A2. This is shown in Fig. 2 (a, 
c, e) where the Fourier transforms of £la(t) are plotted 
for a — {1/3, 2/3, 3/2}27r/A 2 and t g = 4tr respectively. 
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FIG. 2. (color online) Fourier transforms of the control fields 
(a, c, e) and populations (b, d, f) of the ground (blue solid), 
first (red dotted), and second excited (green dash) in a simu- 
lation of a NOT gate with a Gaussian amplitude pulse for a 
d = 5 SNO. In (a, b) a = 27r/3A 2 , (c, d) a = 4tt/3A 2 , and 
(e, f ) a = 371-/A2 and the gate time is taken to be 4<r. 

To quantify this error, we use the gate fidelity averaged 
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over all input states existing in the qubit Hilbert space, 



dip Tr 



t^ideal 



(3.3) 



where S(p) is the actual process in the full Hilbert space. 
Using an argument similar to Ref. [39] gives 



E 



Tr 



j=±x,±y,±z 



UideBlPjU? degl £(pj 



(3.4) 



where pj are the six axial states on the Bloch sphere, 
and ?7idcai is defined in Eq. (3.1). To demonstrate the 
errors arising from Gaussian shaping we consider a d = 5 
SNO and numerically calculate the gate error (1 — F g ) 
for cr = {1/3, 2/3, 3/2}2vr/A 2 and t g = 4a. We hnd gate 
errors of 0.198, 0.0160, and 0.0030 respectively. 

To understand these error values we plot the popula- 
tions of the first three levels in Fig. 2 (b, d, f). The 
ground state populations are given by the blue solid line; 
the red dotted line shows the first excited state; the 
green dashed line is the second excited state. We ob- 
serve that for the shortest gate, o = 27r/3A 2 , the error 
after the pulse is mostly residual population of the third 
and higher level. This is what we refer to as the leakage 
error. For longer gates, e.g. a = 47r/3A 2 , the resid- 
ual population does not account for the calculated error. 
This error is mostly a phase error resulting from the fi- 
nite population of the third level during the pulse. Even 
though the final state is restricted to the computational 
levels, the admixture of the third level leads to a phase 
shift on the second level, resulting in a net phase error at 
the end of the pulse. At the longest time when the pop- 
ulation of the third level is nearly negligible there is still 
a large gate error. From these results we conclude that 
Gaussian shaping is of limited performance even if the 
pulse bandwidth is somewhat smaller than A 2 . Thus a 
more advanced pulse is needed. In Ref. [30] , we provided 
a simple scheme and in the next section we will review 
this and generalize the result. 



IV. SIMPLE ADIABATIC CONTROL PULSES 
A. General procedure for DRAG solutions 

To go beyond the Gaussian control methods presented 
above, we introduce the Derivative Removal by Adia- 
batic Gate (DRAG) technique, generalizing the result of 
Ref. [30]. In this technique we want to find a time- 
dependent unitary transformation A(t) such that the ef- 
fective Hamiltonian 



JTeffW = A\t)H(t)A(t) + iA\t)A{i). 



(4.1) 



has the form H q h(t) © H rest (t) where H q b(t) is a Hamil- 
tonian in the effective qubit subspace and H rest (t) gener- 
ates evolution in the rest of the system. The direct sum 
form implies that if the state of the system starts in the 



qubit subspace, then it remains there during and after 
the pulse. We impose the additional requirement that 
A(0) = A(t g ) = 1, i.e. the frame transformation van- 
ishes at the boundaries. Both conditions together imply 
that we have decoupled the computational subspace from 
the leakage subspace since the qubit subspace is mapped 
back onto itself by the end of the pulse. In Ref. [30] we 
restricted the problem to a d = 3 system and found a 
solution to the problem of generating a NOT gate. It is 
straightforward to generalize this work to design pulses 
that account for higher order corrections, larger dimen- 
sional embeddings, general target gates, and an expanded 
library of control fields by considering more complicated 
frame transformations to that considered in Ref. [30]. 
We now classify valid choices for the transformation A(t) 
that satisfies the above constraints. 

To simplify the following arguments we consider a ba- 
sis for the Lie algebra u(d), as opposed to su(d), with 
elements {a* k ,a y k ,Hi}. Here 1 < I < d — 1 and 
< j < k < a. With respect to this basis we wish 
to implement a Hamiltonian in the qubit subspace of the 
form 

H qh (t) = l - [M*Ki + M*Ki] + \h z {t) (n - n x ) . 

" (4-2) 

Physically this amounts to setting conditions on the con- 
trol fields and the frame transformation such that 

■&[J?eff(tK,i] = h x (t), (4.3) 

Trli^KJ = h y (t), (4.4) 

Tr[H eS (t)(U -U 1 )] = h z (t), (4-5) 

Tr[iZeff(*)CTo,fc] = for 2 < k < d - 1, (4.6) 

Tr[H eS (t)a y o k ] = for 2 < k < d - 1, (4.7) 

Tr[iJ eff (i)of ife ] = for 2 < k < d - 1, (4.8) 

Tr[H cfi (t)a y l k ] = for 2 < k < d - 1. (4.9) 

The three first equations define the controls in the qubit 
subspace whereas the last four impose the condition of 
no leakage. There are no further conditionsso that the 
dynamics strictly inside the leakage subspace remains ar- 
bitrary. 

To find a scheme that satisfies Eqs. (4.3) - (4.9) we 
write A(t) = cxp[— iS(t)] where S(t) is an arbitrary Her- 
mitian operator that we decompose as 

S(t) = s ^(t) u 3 + E 8 *J*V)<>j,k + E a ™.fe(*Kfc> 

j = l j<k j<k 

(4.10) 

and by assuming a power series in a small parameter e, 
we can write each element as 



(4.11) 



n=l 



where a — x 7 y, or z. This ensures that the transforma- 
tion is perturbative with respect to the parameter e. We 
take e = l/i s A as a small parameter, implying that, for 
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the fast gates we are interested in, Oq/A will be of order 
e. Furthermore, since we are interested in unitary oper- 
ations, we will define a dimcnsionless time and Hamilto- 
nian whereby j* 9 H(t)dt = ^ H{t)dt (so H(t) = t g H). 
Doing this allows us to write the Hamiltonian for the 
system in Eq. (2.8) as 



H{t) = -H + 8{t)H z + ^H x + ^H y , (4.12) 



C the following constraints on the frame transformation: 



where 



d-1 

3 = 1 
d-1 

d-l 



(4.13) 
(4.14) 
(4.15) 
(4.16) 



Here S(t), £l x (t), and & y (t) are dimensionless versions of 
the previously defined matrix elements scaled by the rule 
S(t) = t g S(t). Note we have dropped the R superscript 
and assume from now on we are in the rotating frame. 

For the control fields we also write a series expansions 



Q x (i) = £e"0(")(i), 

CO 

n w (t) = 5> n n(">(t), 

11=0 
oo 



(4.17) 
(4.18) 
(4.19) 



n=0 



and by the results presented in appendix C, the con- 
straints Eqs. (4.3) - (4.5) can be rewritten as 

^ n) (*) = h x n) {t) TriH^Ut)^}, (4.20) 
fi(«)(t) = h y n \t) Tr[^ ra (i)<i], (4-21) 
= TriH^Ut) (n - Hi)] - fc<»>(i), (4-22) 



where H^ Ta (t) contains terms generated by the lower 
orders of the transformation and the controls. It is a 
rather complicated expression which can be derived fol- 
lowing the procedure in appendix C and the first few 
orders are listed in Eqs. (CIO) - (C13). For the leakage 
constraints Eqs. (4.6) - (4.9), we also derive in appendix 



,("+!) 



(*) = 



A 2 
2A fe 

A 2 

2A fc 

A 2 



Tr[iJ c ( l a (t)< fe ], 



Tr[^ ra (i)<J, 



(4.23) 
(4.24) 



„(«+!) 



2A« 



A 2 
2A fc 



n^ n) (t)AiJ fc ,2 



(4.25) 



ni n) (t)AiJ fc ,2 



(4.26) 

The coefficients s^^t), s y n ^l\t), and are free 

parameters in our theory. Choosing a different functional 
form for them result in different DRAG solutions. In the 
next section we will give some practical examples for this 
choice. 



B. Zero, first, and second order DRAG solutions 

1. Zero order solution 

For definiteness we choose the target Hamiltonian to 
be Q(;(t)<Toi/2 corresponding to rotations around the x- 
axis. Rotations around the y-axis will follow a similar 
procedure, and z-axis rotations are trivial. For this target 
Hamiltonian we require h x (t) = tgClcit), h y (t) — 0, and 
h z {t) = 0. For simplicity, we take h x °\t) = tgClcit) and 



h ( x n \t) = for n > 0, and wc take h { y n) (t) = h z n '(t) = 
for n > 0. From Eq. (C10), the zeroth order expression 

#extra W i s zer0 - This implies that the control constraints 
in Eqs. (4.20) - (4.22), for order zero, are 

n x °\t)=t g n G (t), fi<°>(t) = o, m(t) = o, (4.27) 

giving control solutions 



n x {t) = n G (t), n„(t) = o, s(t) = o. 



(4.28) 



These are the controls used in Sec. Ill and here we will 
use them as a benchmark for the higher order solutions. 
In Fig. 3, we plot the error, 1 — F g with F g given by 
Eq. (3.4) (blue dashed line), between a NOT gate an a 
unitary from the control field given by Eq. (3.2) with 
A = 7r and t g = 4er for a SNO with d = 5. In this 
figure it is clearly seen that the error associated with 
these controls is quite large; for fast gate times this error 
is unacceptable for quantum information processing, and 
long gate times will have additional error arising from 
decohcrcncc. 



2. First order solution 



To determine the first order solutions we need 
i^extraWi which requires determining the frame transfer- 
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0.5 1.0 1.5 2.0 2.5 
pulse width a [units 27T/A2] 

FIG. 3. (color online) Gate error for the implementation of 
a NOT gate in a d — 5 SNO as a function of a, with t g = 
4<r, and a Gaussian shaped pulse. The blue dotted line is 
the zeroth order solution. The black dash line is the first 
order Z-only correction. The red dash-dot-dot line is the first 
order Y-only correction. The green dash-dot line is the first 
order correction from the controls presented in Ref. [30]. The 
purple solid line is for the optimal first order correction. 



mation conditions for S^(t). From Eqs. (4.23) - (4.26), 
these are 



s S, fe (*) =o, = -\it g n G (t)5 k , 2 /2. 



(4.29) 



Using Eq. (Cll) for i?extra(*)> the first order corrections 
to the control fields are 

^ 1) (*) = 2igJ ll (t), (4.30) 

ni 1 Ht) = 2sl%(t)- s ^(t)t g n G (t), (4.3i) 



(4.32) 



Here we see that there is a continuous family of DRAG 
pulses; however, in this section we will consider four 
particular solutions. In all of these solutions we take 
s^' j (t) — as it has no influence on our choice for 

n£\t) and 5^(t). 

The first solution we consider is one where the control 
field Cl y (t) — 0. This is achieved by setting 4 lW = 

2 s ( yl A (t)/t g n G (t), result! 



mg m 

t g n G (t) tJf G {t) 

+ 2 S « 1 (t)i a n G (i) + ^M), 



(4.33) 



The simplest solution that satisfies (0) = (t g )=0 
is s^o 1 (t) = 0. In this case the controls become 



n x (t) =n G , n y (t) = o, 6(t) 



4Ao ' 



(4.34) 



For the SNO considered in Fig. 3, the error for this con- 
trol set is plotted as the black dashed line. It clearly has 
a much lower error then the standard Gaussian ampli- 
tude modulation control, and we will refer to this as the 
Z-only correction. 

The second control solution we consider is when the 
control field S(t) = 0. This is achieved by setting 

41l(*) = -41iW/ 2 ^G(i) - A?t ff fi G (*)/8, which re- 
sults in 



a)/ 



tgSl G (t) tgto%(t) (4.35) 

- -S(*)*.M*)- 

Again, the simplest solution that satisfies S^(0) — 

S^K^g) = is 4i(*) = 0. In this case, the controls 
become 



n x (t) =n G {t), n y (t) 



4A 2 : 



8{t) =0, (4-36) 



and, for the SNO considered in the numerical demon- 
stration, the error for this control set is plotted in Fig. 3 
as the red dash-dot-dot line. Its error rate is lower then 
both the standard Gaussian controls and the Z-only cor- 
rection. This is the control procedure used in Ref. [32] 
and Ref. [33], where it was referred to as simple DRAG 
and half derivative respectively. Here we will refer to this 
as simply y-only correction. 

The third control solution we consider is what we refer 
to as the optimal first order solution. This is achieved 

( 2) 

by minimizing the elements in H^xtraW such that the 
second order corrections to the control fields are zero. 
From Eq. (C12), we require calculating the second order 
frame transformations. From Eqs. (4.23) - (4.26) these 



,(2) 



1 



S y,0,2 - -2^^ G ^ 1 ^f^ G + S a!o,l)> 



b x,Q,2 ~ 2 9 G 1 z,0,l' 
b y,l,2 — A l b x,0,n 



ot 1 ) 



(4.37) 

(4.38) 
(4.39) 



(4.40) 

Using these expressions and requiring that the matrix 
elements of H^ x ' tla (t) are zero in the qubit subspace (and 
elements coupling to the qubit subspace are zero) results 
M ~ 4!iM = and sW^) = -i ff ft G (t)Ax/4. 
Substituting these into Eq. (4.30) gives the controls fields 



") s ;,o,i(*) 



n x (t) =n G (t), si y {t) = - 

*(*)=^[A?-2A l] . 



2A 2 



(4.41) 
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This optimal first order solution is plotted in Fig. 3 as 
the solid purple line. Its error is substantially lower the 
the other first order correction methods. 

Finally for completeness we also present the first or- 
der DRAG solution presented in Ref. [301. This occurs 



when we choose s^ i(t) = s^\(t) 
— Q<3(t)/2, resulting in 



n x (t) =n G (t), n y (t) 

S(t) 



GV '[X\ -4]. 



order. It is for this reason that we argue this is the best 
solution for implementing a DRAG correcting pulse. 

We have not proven the DRAG solution to be optimal, 
and it seems likely that a better second order solution 
exists. To find the optimal solution, the matrix elements 
= and s^Q 1 (t) = of H^ tIll (t) in the qubit subspace and the elements that 
couple to it must be minimized. Wc have not computed 
this solution due to the complexity of H^lr^t) as well as 
the already outstandingly low error of the original DRAG 
solution. 

(4.42) 



' A 2 ! 



4A, 



This solution can be intuitively derived from an inter- 
action picture; however, there is nothing optimal about 
this choice. In Fig. 3, the green dash-dot line shows how 
the error scales with this control set. We see that for 
the first order solution, the error is larger than both the 
optimal and the Y"-only correction methods. That is, the 
first order solution of Ref. [30] is not optimal. 



3. Second order solution 

The higher order solutions become impractical to solve 
in generality because the number of terms grows quickly 
with increasing order. However, we can easily find the 
second and higher order corrections to the different first 
order solutions. To do this, it is simplest to use a com- 

(2) 

puter algebra system as the expressions for H^ x l ra (t) and 

-^oxtraW are rather involved [see Eqs. (C12) and (C13)]. 
We find that the corrections to the above four cases only 
change the il x (t) field. In the Z-only case, O x (t) becomes 



n x (t) = n G (t) + 



xln G 



8A% 



in the Y-only case, fi x (t) becomes 



n x (t) = n G (t) 



A?(A? 



G 



32A2 



and the control set presented in Ref. [30] gives 



n x {t) = n G (t) + 



(A?-4)ft G 



8A| 



(4.43) 



(4.44) 



(4.45) 



To demonstrate these correction we plot in Fig. 4 the 
second order solutions. The line marking and colors are 
the same as in Fig. 3 with the exception that they now 
refer to second order solutions (all except the blue dot- 
ted line, which remains the zero order solution, and the 
purple solid line, which is the first order optimal solu- 
tion). We see that the second order only makes small 
improvements to the F-only (red dash-dot-dot) and Z- 
only (black dashed) first order solutions. Remarkably, 
the original DRAG scheme from Ref. [30] (green dash- 
dot) is improved substantially when corrected to second 




0.5 1.0 1.5 2.0 2.5 
pulse width a [units 2n/A2] 

FIG. 4. (color online) Gate error for the implementation of a 
NOT gate in a d = 5 SNO as a function of o, with t g = 4<r, 
and a Gaussian shaped pulse. The blue dotted line is the 
zeroth order solution. The black dash line is the second order 
Z-only correction. The red dash-dot-dot line is the second 
order Y-only correction. The green dash-dot line is the second 
order correction from the controls presented in Ref. [30] . The 
purple solid line is for the optimal first order correction. 



C. Numerically optimized first order solutions 

Given the ease of implementing the first order solu- 
tions, in this section we consider the problem of numer- 
ically optimizing a value for the control fields with the 
following ansatz 



n x (t) = a n G (t), n v {t) = -p 



n G (t) 
A 2 : 



5(t) = 7 



(4.46) 



So, 



where a, (3, 7 and So are fit parameters. We consider a 
SNO with A2 = — 2-7T, Xj-\ — y/j, d = 5, and a control 
field given by Eq. (3.2) with A = it and t g = 4cr (same 
as before). In Fig. 5 we plot the gate error as a func- 
tion of a for different optimizations. The optimization 
procedure was done with Mathematica with a working 
precision of 10. In Fig. 5 (a) we consider the case when 
So = 0, and we find that optimizing the weighting of 
the control fields only improves the first order solutions 
slightly. This is expected as the second order solutions 
require different functional forms for the controls. How- 
ever, when we allow So to be non-zero, we find some 



interesting results. For the numerical parameters con- 
sidered, we find that implementing a time varying 8(t) 
(7 7^ 0) does not lead to any improvements. This is seen 
in Fig. 5 (b) where we show that the error arising from an 
optimized Gaussian with added constant detuning (blue 
dotted line) is approximately equal to the optimized Z- 
only correction with an added constant detuning (black 
dashed line). Furthermore, the optimized Y-only correc- 
tion with an added constant detuning (red dash-dot-dot) 
is approximately equal to the optimal first order solution 
with an added constant detuning (solid purple). We also 
find that for the solutions with the derivative for the Y- 
control (solid purple and red dash-dot-dot) the gate error 
is much lower then in the other cases (blue dotted and 
black dashed). This gate error is approximately equal to 
those found with the second order corrections from Fig. 4 
(green dash-dot line). We conjecture from these numerics 
that the optimal DRAG- like solution can be obtained by 
applying a pulse to the x-axis (and its derivative to the y- 
axis) with a frequency that is not equal to the transition 
frequency of qubit. 




pulse width a [units 27U/A2] 

FIG. 5. (color online) Gate error for the implementation of a 
NOT gate in a d = 5 SNO as a function of a, with t g = 4a, and 
a Gaussian shaped pulse. Panel (a) is for optimized first order 
solutions. The blue is for optimized a with /3 — 7 = So — 
(zeroth order solution). The black dashed is for optimized a 
and 7 with /3 = So = (Z-only solution). Red dash-dot-dot is 
optimized a and /3 with 7 = So — (Y-only solution). Purple 
solid line is for optimized a, /3, and 7 with So — (optimal 
first order solution). Panel (b) is the same as (a) but with So 
being optimized. 



V. LEAKAGE EXTENSION 

In this section, we show that the DRAG technique can 
be applied to systems with more then one leakage tran- 
sition. To show this we consider the two cases shown 
in Fig. 6. We note that these are arbitrary examples, 
and the theory is more general then considered here [see 
appenix C]. 



(a) . (b) 

• A 




FIG. 6. (color online) Energy level diagrams for systems that 
can also be model by our theory. In (a) we consider the case 
when the qubit has leakage from both its |0) and |1) level. In 
(b) we consider the case when the qubit has leakage from its 
1 1) level to more then one higher level. 



A. Leakage from both logical states 

The first case we consider is a qubit defined in the in- 
termediate states of an anharmonic oscillator [Fig. 6 (a)]. 
This situation is important if the anharmonic oscillator 
is going to be used for qudit logic, as done in Ref. [40], or 
for state tomography of the qudit, as done in Ref. [41] . In 
this case, we rewrite the free and coupling Hamiltonian 
as 

N 

H h .= (ju> + A j )Il j , (5.1) 

j=-N 

N 

H ct (t)=E{t) Vi^-l, (5-2) 

j=-N+l 

where N = (d — l)/2 and again we take Ao = Ai = 0. 
Moving to a interaction frame similar to Eq. (2.4) (the 
sum range is change to be consistent with the above) we 
find a dimensionless rotating frame Hamiltonian equiva- 
lent to Eq. (4.12) with 

N A 

H = J2 ( 5 - 3 ) 

j=-N >7 iO,l 2 
N 

j=-N 
N 

j=-N+l 
N 

Hy= £ 'V. i< (5.6) 

3=-N+l 

Using the results of appendix C, the zero order dimen- 
sionless controls are the same as Eq. (4.27), implying 
that the zeroth order controls are given by Eq. (4.28). To 
find the first order corrections we follow a similar proce- 
dure to Sec. IV. The frame constraints from Eqs. (4.23) 
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and (4.24) become 



„(»+!) u\ _ _ A 2 



■\,,i./. l" = -^Tr[F^ ra (tK i)t ]- 



2A_i 



(5.7) 



(n+i; 



(*) 



A2 

2Az 



lt[^ ra (t)< fc ] 



A 2 n< n) (i)A_i4,- 



2A_ 



(5.8) 



while Eqs. (4.25) and (4.26) remain the same. With 

-HextraW gi ven by Eq. (Cll) we find the first order frame 
transformation to be 



L(*)=o, »SS.*(*) = - 



,<i) 



U(*)=°. »£U*) = - 



A_ 1 A 2 ^^ G (t)^,_ 
2A_i 

Alt g fic?(i)£fc,2 



and the dimensionless first order control fields are 
n«(t) = 2iW 1 (i), 

^ 1) w = 24ii(*)-4!l(^ G (*), 
5 (1) (*) = 4!i(*)+241o,i(*)*^w 

''/ :i gW/^ A2 -A 2 



■4 l A i~ AT7 A -i 



(5.9) 



(5.10) 
(5.11) 
(5.12) 

(5.13) 



From this we find the control fields for Z-only correction 



shown as a function of a for the same Gaussian shaped 
pulse as considered in Sec. III. Here we observe that 
the DRAG technique improves the gate fidelities sub- 
stantially when compared to the zeroth order solution 
and has error rates comparable to that of the case with 
only one leakage channel. 




0.5 1.0 1.5 2.0 2.5 
pulse width a [units 271/A2] 

FIG. 7. (color online) Gate error for the implementation of a 
NOT gate when there is leakage above and below the qubit 
subspace. The system considered is explained in the text. 
The blue dotted line is the zeroth order solution. The black 
dash line is the first order Z-only correction. The red dash- 
dot-dot line is the first order F-only correction. The purple 
solid line is for the optimal first order correction. 



n x (t) =n G (t), n„(t) - 0, 



S(t) 

the F-only are 

n x (t) =n G (t), n y (t) 
S(t) =0, 



(5.14) 



n G (t) 



a: 



(5.15) 



and the optimal first order control field corrections (after 
minimizing H^ ia (t)) are 

(l G (t) 



S(t) 



2A 2 
4A 2 



i2_i_ A 1 A 2 



A? 



A_ 



(5.16) 
(5.17) 

(5.18) 



To numerically demonstrate an improvement over the ze- 
roth order solution we consider a SNO (of d = 6) where 
we want to control the 2 — > 3 transition. In this case, we 
relabel j = 2 to and so on, rescaling the coupling so 
that the new — > 1 transition is unity. This results in 
setting A 3 = 3A 2 , A_i = A 2 , A_ 2 = 3A 2 , and A = 1, 
A x = ^4/3, A 2 = y^5/3 A_r = A_ 2 = y/l/3. In 

Fig. 7, the gate error for implementing a NOT gate is 



B. Leakage from the excited state to more then 
one auxiliary level 



In the second case, shown in Fig. 6 (b), we consider a 
qubit made from the lowest 2 levels of a system which is 
coupled to many other transitions, all transitions having 
only a small energy cost (approximately A 2 ). This is an 
interesting example as it shows how this theory can be 
easily generalized. In this case we rewrite the free and 
coupling Hamiltonians as 



(5.19) 



H h = will + ^(2w + Aj)TLj, 
H ct (t) = £{t) (a <i + £ VKi) ■ (5-20) 



3 =2 



To eliminate the fast degrees of freedom we move to a 
rotating frame and make the standard rotating wave ap- 
proximation. The procedure is similar to Sec. II B with 
the replacement of Eq. (2.4) by 



d-l 



R{t) = exp (— iuj d t) III + ^2 exp (~ i2uj dt) lij- (5.21) 



J=2 
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This results in a dimensionless rotating frame Hamilto- 
nian equivalent to Eq. (4.12) with 



d-l 



Ha 



3=2 



i=2 



d-l 



d-l 



3=2 3=2 

Again, we find the zeroth order controls given by Eq. 
(4.28). To find the first order corrections, we follow a 
similar procedure as in Sec. IV, finding that the frame 
constraints Eqs. (4.23) and (4.24) remain the same, but 
Eqs. (4.25) and (4.26) are changed to 



2A fc 



-£S* 1) (*) = ^[flSL(*KJ 



2A fc 

(5.22) 

^ n) (t)A fc -iA 2 
2A fc 

(5.23) 

From the above with i?^ ra (t) given by Eq. (Cll), we 
find the first order frame transition to be 

*£o,*(*)=0, »$,*(*) =°> 

S (1) (t)-0 S (1) (t)~ A ^Wt)A 2 (5-24) 
This gives the dimensionless first order control fields 

=2A$, 1 (t) ) 

=2AS, 1 (*)-«S(t)* fl nG(*) J 

« (1) (t)=ii3(t)+4,i(*MG(*) (5 - 25) 

From this, we find the control fields for the Z-only cor- 
rection are 



n x (t) =n G (t), n y (t) = o, 



S(t) 



E 



fc=2 A fc 



(5.26) 



for the Y"-only correction are 

n a (t) =n G (t), ^W = -^Et2%: 
*(*) =o, 



(5.27) 



and the optimal first order control field corrections (after 

(2) 

minimizing i?extra(*)) are 



<>,(/>=<><;(/>• - n ^ A fci 



1 "2"fc-l 

2 A? • 



5(t) 



' 4A, 



^fe=2 A? V Z^fe=2 A? 



(5.28) 



We note that these solutions are identical to the previous 
solutions with a single leakage channel where Ai = A, 



d _i A|A|_i 
2^fc=2 A? ' 



(5.29) 



To numerically demonstrate an improvement over the 
zeroth order solution, we consider the implementation of 
a NOT gate for a d = 6 system with Aj = 1 for all j and 
A3 = 2A 2 , A4 = 3A2, A5 = 4A2 (note these are different 
parameters from the anharmonic oscillator considered in 
Sec. IV). The results are plotted in Fig. 8, where again, 
it is clearly seen that the DRAG technique improves the 
zeroth order solution. 




0.5 1.0 1.5 2.0 2.5 
pulse width a [units 271/A2] 

FIG. 8. (color online) Gate error for the implementation of 
a NOT gate when there is many leakage transitions for the 
excited state. The system considered is explained in the text. 
The blue dotted line is the zeroth order solution. The black 
dash line is the first order Z-only correction. The red dash- 
dot-dot line is the first order K-only correction. The purple 
solid line is for the optimal first order correction. 



VI. CONCLUSIONS 

In this paper we have presented a general technique 
for designing simple controls fields for single qubit uni- 
tary operations in weakly non-linear oscillators, which 
we refer to as the DRAG (Derivative Removal by Adia- 
batic Gate) technique. We first consider a qubit formed 
by the two lowest levels of an anharmonic oscillator with 
only the one photon transition elements being non-zero. 
In this system, the largest source of error for fast gates 
(small gate times) is leakage from the |1) state to the \2) 
state. Our technique provides a simple control method- 
ology that perturbatively removes this error, thereby al- 
lowing high fidelity single qubit gates. The essential idea 
of this method is to apply a y-field that is proportional 
to the derivative of the original control pulse. It contains 
the DRAG solution presented in Ref. [30] as well as a 
large collection of control pulses that also correct for this 
leakage error, which, for example, require fewer control 
fields. The lowest error obtained requires this derivative 
correction as well as a frequency shift. Furthermore, we 
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show that this methodology can be easily extended to 
other weakly non-linear oscillators with more then one 
leakage transition. 

In this paper, we considered Gaussian pulses due to 
their favorable spectral properties, however, our theory 
is independent of the initial form of the pulse shape. A fu- 
ture research direction could be to find the optimal pulse 
shape for the leakage problem. Furthermore, while the 
pulses found here are robust to first order in variations in 
the control parameters it would be interesting to search 
for pulses with higher order robustness properties. 
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Appendix A: Typical functional forms for \j 

In this appendix, we discuss the different possible val- 
ues Xj for the realistic models of direct driving [17] and 
the cavity [29, 42] or circuit [43, 44] QED architecture. In 
cases where the system is a SNO and is driven directly, Xj 
is well approximated by the harmonic oscillator matrix 
elements, namely Xj w y/j. This, for example, occurs 
for the phase qubit when it is driven by a time-varying 
bias current [17] and the transmon when it is driven by 
time- varying gate voltage [45]. Essentially, these systems 
are very nearly harmonic oscillators and the controls are 
almost proportional to the quadrature operator. 

In the cavity or circuit QED architecture, the anhar- 
monic oscillator is coupled to a resonator and is control 
by driving the resonator far off-resonance. In this case, 
Xj can take on essentially any value. To see this, wc 
start by writing the full Hamiltonian for the multi-level 
anharmonic oscillator and resonator as 



Hj C =Lo r a)a + EC? W + Aj)Uj 



3 = 1 



d-l 



(Al) 



+ E»-^(lJ- 1 >0V + li>0'-i|o) 

i=i 



again Ao = Ai = 0. This is a generalized Jaynes- 
Cummings Hamiltonian [46] with u r being the resonator 
frequency, and gj^ being the vacuum Rabi coupling for 
the j to k levels. Again we have assumed that the an- 
harmonic oscillator only allows one photon transitions. 

If j — w r | ^> \gj-i,j\ for all j where u'j-ij = 
u> + Aj — Aj_i then diagonalization of Eq. (Al) can be 
performed to lowest order in gj_ij/(uij-ij — uo r ) by the 
canonical transformation Hf c = HjqD where 



D 



cxp 



d-l 



E , gj ~ 1,j (a f |j-l)01-4?)(j-l| 

— d J ■ - ■ — CU- 



3 = 1 ^-h3 



to give 



d-l 



[LO r 



Xoi)a f a + C? w + A j + Xj-i,j) Hj 



3 = 1 



(A2) 



(A3) 



E 5 ^-' + *j-u - xjj+i)° tan J- 

3 = 1 



Here Xj-i,j = 



3-1,3 



is the Lamb shift in- 



duced on the anharmonic oscillator by the resonator and 
the last term in Eq. (A3) is the ac-Stark shift [47]. As- 
suming that the dressed cavity is in vacuum then Eq. 
(A3) is well approximated by Eq. (2.1) with 



A, 



to -- 



u + Xo.i 

= A 3 + Xj-1,3 



3X0,1 



(A4) 
(A5) 



and the tilde implies dressed values for the transition 
frequency and anharmonicity from the resonator. 

As stated above, the control is usually through the 
resonator and is represented by the Hamiltonian 

H dl (t) = e(t)(a + a t ) (A6) 
which under the transformation Eq. (A2) becomes 



d-l 



ff£(t)=e(t) a + at+E-7 



9 3 -1,3 



j=l 3-1,3 r 



'3-1,3 



(A7) 



Assuming that e(t) is a sinusoidal with a frequency close 
to the qubit then Eq. (A7) is well approximated by Eq. 
(2.2) with 



£{t) = 



9o,i 



W 1 — u r 



<t) 



93-1,3^0,1 ~ u r) 
90,1 K-l,j - Ur)' 



(A8) 
(A9) 



To demonstrate the functional form of A l7 we will 
assume the anharmonic oscillator is in the SNO limit: 



o'j_ hj = u + (j - 1)A 2 , and gj_ hj 
~!q. (A9) becomes 

VJ 



case, Eq. (A9) becomes 
Aj_i = 



\ + {j-i)A 2 /{u-uj r y 



Vj9o,i- In this 



(A10) 
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and is plotted in In Fig. 9 as a function of the bare 
anharmonicity in units of (u — uj r ). Here, it is clearly 
seen that Ai is of not equal to \[2 and actually changes 
sign at the point (w — uj r — — A2). This point inval- 
idates the diagonalization for the second level and can 
not be treated under this model. Away from this point 
the values obtained from this model approximate the real 
situation. This was confirmed in Ref. [32] where the ex- 
perimental value for Ai for the operation point used was 
found to agree with Eq. (A9). However, we propose that 
Xj should be used as a fitting parameter in any exper- 
iment as effects such as higher modes of the resonator 
and higher order perturbation will result in additional 
corrections to this value [48]. 




[.0 -0.5 0.0 0.5 1.0 
A2/(C0-(D r ) 



FIG. 9. (color online) Ratio of coupling strength of the |1) — >• 
1 2} transition to the |0) — > |1) transition as a function of 
anharmonicity for control through a resonator (solid blue) 
and direct drive (dashed red). 



and by using Eq. (2.5) we find an effective Hamiltonian 
of the form 



d-l 



c -iu> A t-i f*6(B)dB ! h|| 



d-l 

x Aj-i - i| e *"**+</o*W«'" +h.c. 

(B3) 

where ft(t) — fl x (t)+iCly{t). Now provided we can make 
the rotating wave approximation this simplifies to Eq. 
(2.8). 

Note Eq. (Bl) may appear that is requires more con- 
trol then is available but this is not true. Using simple 
trigonometric identities, this can be rewritten as 

£ (t ) =n' x (t ) cos [u d t] + Of (t) sin [u d t] (B4) 



where 



Q' x (t) = tt x (t) cos 

n' y (t) = Cly(t) COS 



S(s)ds 



S(s)ds 



+ £l y (t) sin 



fl x (t) sin 



5(s)ds 



(B5) 



S(s)ds 



(B6) 



This is known as phase ramping [49]. 



Appendix C: General leakage removal 



Appendix B: Detuning control for fixed frequency 
qubits 

As a starting point for our control methods, we argue 
that the system can be explained by a Hamiltonian of the 
form displayed in Eq. (2.8), with 5(t), £l x (t), and £l y (t) 
are independent control parameters. We argued that 5(t) 
is achieved by tuning the qubit frequency. Here, we show 
that it can also be achieved for fixed frequency qubits 
by changing the control field. In place of Eq. (Bl), we 
assume that input field is of the form 



£{t) =Q x (t)cos 



ujdt + / S(s)ds 


t 



fly (t) sin 



uj d t + / S(s)ds 




(Bl) 



Now following the same procedure as in Sec. II B with 
£{t) replaced by Eq. (Bl) we can define a frame trans- 
formation R(t) by 



R(t) 



d-l 

J^exp 



UJ d t 



5{s)ds 



n, 



(B2) 



In this appendix, we give the general theory for finding 
adiabatic pulses for controlling subsystems of a larger 
system. The systems that will obey this theory are those 
in which transitions out of the subspace occur an energy 
cost A. Furthermore, we assume this cost is the largest 
rate in the problem so that we can form a perturbation 
expansion around e = l/t g A, where t g is the gate time. 
To be more specific, imagine a Hilbert space that can be 
partitioned into two subspaces, H = ^control © ^leakage, 
where we wish to perform some quantum gate solely in 
the control subspace, but with control Hamiltonians H m 
that act on the complete Hilbert space. That is, the 
general Hamiltonian for this system is 



H(t) = H + J2u m (t)H m , 

m 

where we require that Hq takes the form 

ff = o® A^^n fe . 



(Cl) 



(C2) 



Here = \k){k\ and we note that H only has support 
on the leakage subspace. Since, A is assumed to be large 
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then Hq is the dominant term in the Hamiltonian, and 
the leakage subspace is off-resonant. It is less clear what 
constraints must be placed on the control Hamiltonians, 
H m . Certainly, if we wish to perform a general gate on 
the control subspace, then they must be controllable on 
that subspace, however, in order to completely eliminate 
leakage errors there are additional constraints, which we 
now determine. 

The general protocol to determine these extra con- 
straints arising from leakage is to first determine the 
Hamiltonian in the control subspace that implements the 
desired gate. That is, find the functions hi(t) of the 
Hamiltonian H sp 

H sp = Ht)Bi (C3) 

where Bi is some orthogonal basis for operators in the 
control subspace of dimension d c , that implement the 
desired gate. Next we define an effective Hamiltonian 
H e ff(t) in the complete space by Eq. (4.1) with the 
only requirement on A being that A(0) = A(t g ) = 1 
so that the frame transformation vanishes at the bound- 
aries. The constraints on this effective Hamiltonian to 
be satisfied are 

Tl[H eS {t)B l ] = h l {t), (C4) 

{j\H eS (t)\k) = for < j < d c and d c < k < d, (C5) 

Here Eq. (C4) constrains the H c s(t) to be equivalent 
to Eq. (C3) and Eq. (C5) cancel the coupling between 
the control subspace and the leakage subspace. If the- 
ses constraints can be satisfied we have found a control 
set {u m (t)} that implements the desired gate with out 
leakage. In practice, this will not be possible but using 



the following procedure we can develop a perturbative 
approach. 

The adiabatic transformation A can be written as 
A(t) = e~ lS (*) and we can write both H(t) and the frame 
transformation S(t) as power series in the small param- 
eter e, as 

oo 

Hit) = -H + J2^H (n) (t). (C6) 

n=0 

and 

oo 

s(t) =j2 s{n) (ty- (C7) 

n=l 

Note the power expansion starts at n = 1 as wc restrict 
the transformation to be a perturbation in e. Rewriting 
Eq. (4.1) in terms of 5, as opposed to A, we obtain 

oo 

H eS (t)=J2(^[iS(t)}) j H(t)/j\ 

+4 (l>mt)/^ (eh)^)/^ , 

(C8) 

where ad [^4] is the linear superoperator defined by 
ad[A]_B = [A,B] and we have used a Baker- Campbell- 
Hausdorff (BCH) and Taylor expansion for the first and 
last term respectively. From this form of H e g(t) we can 
isolate powers of e and arrive at 

H$(t) = H£Ut)+H^(t)+*lS {n+1) (t),H ] (C9) 

for n > 0. Here H^ ra (t) is a rather complicated expres- 
sion which for the first few orders is 



H^L(t) = (CIO) 
H^(t) = i[SW(t),H^(t)] - [S^(t),[S^(t),H ]}/2-S^(t). (Cll) 
lf£L(t) =i[SV(t),HW(t)] +i[SW(t),HW(t)] - [S^(t),[S^(t),H^(t)]]/2-[S^(t),[S^(t),H }}/2 

-[S< 2 >(i), [SW(t),H ]]/2 - i[SW(t), [SW(t), [S^(t),H }}}/6 + t[S^(t),S^(t)}/2 - S^(t). (C12) 
H%LV) = i[S^\t),H^(t)}+i[S^(t),HW(t)}+i[S^(t),H^(t)} [SW(t),[SW(t),HW(t)]]/2 

-[SW(t), [S^(t),H(°\t)]}/2 - [SM(t), [S^(t),H^(t)}]/2 - i[SW(t), [SW{t), [S^(t),H^(t)}}/6 
-[SW(t), [S^(t),H }}/2 - [S< 3 >(t), [S^(t),H }}/2 [S< 2 >(t), [S^(t),H }}/2 

-i[SW(t), [SW(t), [S^(t),H }}}/6 i[S< 2 >(t), [sW(t), [S^(t),H a }}}/6 i[S^(t), [S™(t), [S^(t),H ]}}/6 
+[S< 1 >(i)JS< 1 >(t)JS( 1 >(t)JS^ 

-[S {1) (t), [S (1) (t), S {1) (t)}/6 - S {3) (t). (C13) 



This expansion bears some formal similarity to the 
Schrieffer- Wolf/van- Vleck/adiabatic elimination expan- 



I 

sion [50] which means that a compact closed expression to 
any order is unlikely to exist. This expression for (t) 
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is useful as it allows us to separate the free variables 
if(")(i) and S (n+1 \t) for the order n expression from 
those that are used to satisfy the conditions in Eqs. (C4) 
and (C5) for order n—l, which only occur in the expres- 
sion for ffgxtraW- To see this more clearly we can used 
the above to rewrite Eq. (C4) for each order n as 



(n) 



Tr[i4ta(^]. 



(C14) 



for < I < d\ where is the n th order expression for 
hi. Furthermore Eq. (C4) for each order n becomes 



A 



= **-01 



\k) (C15) 



for < j < d c and d c < k < d. Here we have used that 



H a has no weight in the control subspace. Since i?extra 
only depends on sjj, tf , and ifC™" 1 ) then Eqs. (C14) 

and (C15) are well formed. The components of S'( n+1 ) (t) 
that haven't yet been defined are free parameters in our 
theory. Choosing a different functional form for the re- 
maining matrix elements gives us the different forms of 
a DRAG solution. It is necessary that the frame trans- 
formation vanish at t — and t = t g , which may restrict 
some of the freedoms we have in assigning value to the 
otherwise unconstrained elements of S^ n+1 \t). At each 
order of approximation, specifying the values of (t) 

and S^ l+1 \t) fully determines H^\t) and so we can 
iterate this procedure indefinitely, though at some point 
the derivatives of S that appear in H c g(t) will likely 
become large, causing the adiabatic expansion to break 
down. 
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